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ANALYSIS OF ELECTROMAGNETIC INTERFERENCE FROM 
POWER SYSTEM PROCESSING AND TRANSMISSION 
COMPONENTS FOR SPACE STATION FREEDOM 


1 Introduction 

The goal of this research project was to analyze the potential effects of 
electromagnetic interference (EMI) originating from power system processing and 
transmission components for Space Station Freedom. The approach consists of four 
steps: 


1. Develop analytical tools (models and computer programs). 

2. Conduct parameterization (what if?) studies. 

3. Predict the global space station EMI environment. 

4. Provide a basis for modification of EMI standards. 

This report along with previous interim reports constitutes the final report for 
work performed during the period January 1, 1990 - June 30, 1992. 

The Device Models portion of the investigation centers on the development of 
computer models which can be used to predict the local electromagnetic fields for 
various power system components or devices, including various transmission line 
configurations. This work was summarized in previous interim reports. More 
recently, there has been interest in the vxB voltage drop developed along the length 
of the space station structure as a result of the motion at velocity v through the 
earth’s magnetic B field. This voltage difference along the structure results in a 
current in the conducting plasma surrounding the space station. An analysis of the 
v x B effect, including plots of the current density distribution for a cylindrical model 
of the space station, are summarized in Section 2 of this report. 

The Radiation Sources portion of the investigation centers on the development 
of computer models to simulate the propagation of electromagnetic waves in the 
ionospheric plasma and the radiation from antennas simulating EMI sources on the 
space station. This work is summarized in Section 3 of this report. 


1 



2. Device Models 


1 Introduction 


An electrodynamic system of the space station freedom in the earth’s geomagnetic field in 
the ionosphere is shown in Fig.2.1. A cylinder type space station is moving at a speed, V ,, 
on its orbit surrounding the earth in the geomagnetic field, B g . The motion of the space 
station through the geomagnetic field results in an electromotive force (emf), V , x B g , along 
its length and causes potential drop across the space station. Because the plasma ionosphere 
environment surrounding the space station is conductive, a current between the space station 
and the plasma is induced. This work presents the governing differential equations of this 
V , x B g problem, the formulation of the 3D finite element (3D-FE) model for this problem, 
with the electric scalar potential (ESP) as the field variable, and the consequent computation 
of the induced current distribution. Due to the infinite boundary of this problem, the 
absorbing boundary condition (ABC) method is mentioned as one possible future work, for 
dealing with this type of open boundary problems. 


2 Derivation of Governing Equations 


For the problem described in the previous section, the total region can be divided into two 
subregions, and (See Fig. 2. 2). Region, fti, represents the structure of the space 
station and is moving at a speed, V 4 , in the geomagnetic field, B g . In region, Q 2 , the 
conductivity <72 equivalent to the plasma environment may be anisotropic. Surface, 5, is 
the interface boundary between region, Qj, and region, on which some electric field and 
current density continuity conditions should be satisfied. Since the potentials will tend to 
zero while the boundary, T, is extending to infinity, the asymptotic (or absorbing) boundary 
conditions (ABCs) can be applied on the outer boundary, I\ 

In region, jji ; 

Because of the V x B component, there are two kinds of electric fields in the region 
fli: an induced electric field, E{, and an external electric field, E e . The induced field, Ei, 
can be expressed as follows: 

E x = V,xB g (1) 

The external field, E e , is a conservative field, so that it can be defined as follows: 

E e = — Vfc (2) 
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where, <f > i is an electric scalar potential (ESP) function in region, Sli . The total resultant 
electric field, E\, is the superposition of these two fields, that is: 


E-i = E t + E e - V 3 x B g - Vfa (3) 

The current density, J x , is proportional to the total electric field intensity, E i, hence 

Ji = aiEi = &i(Ei + E e ) = — - (V, x B g )) (4) 

where, < 7 i is the conductivity of region, Oj, which is homogeneous. The following current 
continuity condition must be satisfied anywhere: 

V ■ Ji = 0 (5) 

Substituting (4) into (5), one obtains: 

V • [a x (V^ - (V, x B g ))} = 0 (6) 

Because both a x and (1/, x B g ) are constant, V • [cr\(V s x B g )] = 0, equation (6) can be 
rewritten as follows: 

v • [<Ti(V^)] = 0 inQ, 1 (7) 

In region, f) 2 : 

Only an external field exists. It can be defined as follows: 

E 2 = V <p2 (8) 

where, <j) 2 is an electric scalar potential function in region, fi 2 - The current density, J 2 , in 
region, fl 2 , can be written as follows: 

J 2 — <y 2 E 2 (9) 

where, <^ 2 is the conductivity matrix of region, fi 2 , resulting from the possible anisotropic 
properties of this region. Using the current continuity condition again, one obtains: 

V-J 2 = 0 (10) 
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Substituting (9) into (10) gives: 


V • \a 2 E 2 ] = 0 


( 11 ) 


Meanwhile, substituting (8) into (11) yields the following: 


V • \p 2 {^74>2)] = 0 inQ 2 


( 12 ) 


On boundary, S : 

From electromagnetic field theory, on the boundary, S, the tangential component of the 
electric field should be continuous and the normal component of the current density should 
also be continuous. These two continuity conditions can be written as follows: 


f E ■ dl = 0 -> (E 2 - E x ) t = 0-> [-V^ 2 + Vfa -F,x B g ] t = 0 (13) 

£ J • ds — 0 — + ( J 2 — J i) n — 0 —> [—cr 2 V<f) 2 — <Ti(— + V 3 x i? 3 )] n = 0 (14) 


where, the subscripts, n and t denote the normal direction and tangential direction on 
boundary, S, respectively. 

On boundary, T : 

An m-th order asymptotic boundary condition [1] can be assumed to be satisfied on 
the boundary, T. This can be expressed as follows: 


B m <f> 2 — 0 


(15) 


where, B m is the m — th order differential operator. 

Governing Equations: 

In summary, the governing equations for the (V s xB g ) problem with boundary coupling 
conditions and asymptotic boundary conditions are as follows: 

V.MV^)] = 0 mfij (16) 
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V • [f 2 (V^)] = 0 


(17) 


in 0,2 


(V^a - V<f) 2 ) t = (V, x B g )t on S 


(18) 

MV&) — 0"2(V<^2)]n = 01 (Vi Bg)n 

on S 

(19) 


and 


B m <p 2 = 0 on T 


( 20 ) 


3 Finite Element Formulation 


Using the Galerkin- Weighted Residual (GWR) method [2], one can obtain the integral equa- 
tions equivalent to the differential equations which were derived in the previous section, 
equations (16) through (20). The integral equations can be discretized by finite element 
methods [2]. 


3.1 The Galerkin- Weighted Residual (GWR) method 

Using the GWR method, the integral equations can be obtained from the partial differential 
equations (16) and (17) as follows: 



WV • (<ri(Vfc))<Kl = o 

(21) 

wv • (f 2 (vfa))dn = o 

(22) 


where, W is an arbitrary scalar weight function. 

From vector calculus, one can write a vector identity as follows: 

V • (f~F) = /V • F + V/ • F =► /V • F = V • (fF) - V/ -F (23) 
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where, / is any scalar function and F is any vector function. Let F = and / = W, 

one can obtain the following useful transformation from (23): 

WV • (aV<f>) = V • (W{oV<j>)) - (VW) ■ {aV<f>) (24) 


Using the relationship shown in (24), one can rewrite (21) and (22) as follows: 

[ wv ■{a 1 (v<t> 1 ))<m= f v • (Wfavfaydn - [ (vw)-{(r 1 v<t> 1 )<m = o ( 25 ) 

J n l J fii Jfti 

f WV • = / V • ( W(a 2 V<t> 2 ))<m - [ (VW) • (^ 2 V^>)dft = 0 (26) 

Jn 2 J&2 


From vector calculus , Gauss 1 theorem gives: 


/ V * Fdtt = i F-ds 
J o Js 


(27) 


where, is a space volume, the closed surface of which is S. The normal direction s is out 
of the volume. 

Applying Gauss’ theorem, equation (27), to (25) and (26), one can obtain the following: 


- / (VW) • (ffiWi )dn + i • ds = 0 

J n i j s i 


(28) 


W (cr 2 ^ ^ 2 ) * ds + $> W{a 2 V<f> 2 ) • ds — f (VVF) • (<72V<^2)dfi — 0 
/s 2 «/r J ^2 


(29) 


where, S x = -5 2 . The second term of the left side in (29) can be dealt with the asymptotic 
boundary condition (ABC) method which will be discussed in future work and here assumed 
to be zero (the nature boundary condition). That is: 


j W(a 2 Vfa) -dS = 0 


(30) 


Therefore one can merge equations (28) and (29) into the following: 


f (VW)-{a 1 V<t> 1 )dn+ l (VVU) ■ (a 2 V<fi 2 )dn = l WMVfc)-a 2 (Vfc)] n efe (31) 
J n t Jn 3 Js i 
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Substituting the boundary condition (19) into (31) yields: 

/ (VW)-( <7iV<^)dft+ / {VW)-(5 2 V<f>2)dSl= <f Wai(V, x B g ) n ds (32) 

yn! ^02 ■ ,5 i 

To simplifying the problem, one introduces some definitions as follows: 


£|fii = 


£|n 2 = °2; 


(33) 


and 


<f> Ini = <f>u <^ln 2 = <^2 

Using the above definitions, one can rewrite equation (32) as follows: 
/ (W) • = l Wa(V, x B g ) n ds 

./n=fiiun 2 Jsi 


(34) 


(35) 


3.2 Discretization by FE Method 

Define an interpolation function of the scalar potential in a finite element which has m nodes 
as follows: 


r = £ JWi 


j=h 


(36) 


where, is a shape function, and <j>, is the potential value at node j. The node numbers of 
a given element are l\, • • • , l m - 

Let us define the weight function in (35) as follows: 

W = Ni i = ( 37 ) 


Applying (36) and (37) to (35), one can obtain the element integral function as follows: 



VNj<t>j)d ' n = £ € Nid(V $ 
% — 


x Bg) n ds 

J * ' ' 


(38) 
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Moving the summation out of the integral, one can rewrite the previous equation as follows: 


f m 0 

El/ 

~t Jn « 


= i Ni<L{V s 
Js* 

i,= h, 


x 5g)n<is 

^2i ‘ ' i /m 


(39) 


This relationship can be written in the matrix form as follows: 


& 1 U 12 ••• 


’ fai 


'K ' 

. Sl mi l 2 • ■ • J 


. fam . 


.bu. 


Hence, the elemental FE equation can be expressed as follows: 


(40) 


S' ■£' = b' (41) 

where, 

s %j = [ (VNi£LVNj)dn 
Ja* 

bi= f NiZ.{V, xB g ) n ds ( 42 ) 

^ 5 J — ^1 y ^2y * * * > 

and 


fa = Au ] 1 ( 43 ) 

The tangential continuity condition on the interface, S, (Vfa — V<f > 2 ) t = (F, x B g ) t , 
can be considered as a potential jump distribution by which the previous FE formulations 
should be modified. In the next section, the potential jump distribution will be discussed. 


3.3 FE Formulations Modified by Potential Jump Distribution 


On surface, S, we define the potentials as follows: 

<f>\s - fa] (fa - <£i}|s = T (44) 

Consider an element, e, which belongs to region, fii, and on interface, S. The total m 
nodes on element, e, are l\,l 2 , • • • , / m -n, ^21" ’ ■> * n which l lt are on the interface, 
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S e . Meanwhile, the interpolation function in such kind of elements can be expressed as 
follows: 


# = if 

>=< 


(45) 


Because of the fact that <f>i = <f> 2 - T on nodes / x , / 2 , • • • , l n > the previous equation can be 
expressed as follows: 

« = Ni* U + E N,(h, -T,) = { E N,h, + E N ’M - E Wi («) 

J=h j-l[ 3-h j=l[ ]=l[ 

Applying the definition of the potentials in (34) and (44), one can rewrite the previous 
equation as follows: 


f=E N Ai - E N > T i 

j=ii ,.i; 


(47) 


w = E vn,*, - E VN > T , 

3-h 


(48) 




For an element which belongs to region, Q x , and on surface, S, one can obtain the elemental 
equation as follows: 


£[/ (VN l2 VN J )dQ}4> ] -t t [f = <f N t a(V 

} =h Jn ' j=i[ Ja l s ‘ 


s X Bg) n dS 


(49) 


i — l \ J ^2 j * * * ? 


E {J a (VN, SL vN j )in)4 > , = jf _ Marty, x b,),<u + E [j n ^Na.vN,)ia}T, , 

j=h 1 j=l[ 1 

® == ^1) ^2) * " ' i 

Rewriting the modified FE elemental equations in matrix form one obtains the following: 
■SJi.il •S/ l ,i2 ■Si 1 .im 





1 

^ ... 

1 


i 

... o-< 

i 


. 4>i m . 


» 

£ 

... ~cT 

l 


(51) 
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where, 


Si } - 

bi = 


: / (VNiS.VN})dn 

J n* , 

<£ N t a(V , x B g ) n ds + ^ s l} Tj 
fSt J=l[ 


J — ^1 j ^2 j * * * ) 

Z = ? 1 , ^2 ? ’ * ‘ j J = 


( 52 ) 


The computation of the potential jump distribution, T, will be presented in the next section. 


4 Computation of the Potential Jump Distribution 

4.1 Surface Element Analysis 

From equation (18), one can see that the potential jump distribution on surface, 5, is 
governed by the following differential equation: 

VT = -(V, x B g ) t = -E vt (53) 

where, T = (<f> 2 - ^Ols 

One can choose a cost functional to be minimized on boundary, 5, as follows: 


r - NE ‘ r - 

F(T) = J s | VT + E vt \ 2 ds = £ j se |VT + Eyt | 


2 ds 


(54) 


where, N E 3 is the total number of the surface elements, and S e is the area of a given surface 
element. Here we assume that for every element in region, fix, no more than one element 
surface is on the surface, S. In this work, first order tetrahedral finite elements are chosen 
and only nodes, l' ,rn and n , can be on surface, S, (see Fig.2.3.). 

Define an interpolation function on surface S e as follows: 


T = + N m ’T m < + N n -T n ' (55) 

where, N t > , N m > and N n > are shape functions on S e , and ,T m ' and T n > are the potential 
jump distribution on the nodes, l , m and n . 

A local coordinate system (u,v, w) should be adopted (see Fig.2.3.). 


T = T(u,v) = iV)'(u,u)Tj» + N m '(u,v)T m > + N n >(u,v)T n 


(56) 
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Here 


"fcr, + d -^T . + • 

8u 1 flu m du “ 


§T_dN LT .dN^r . W.- T 

8v _ a» T '' + 9v Tm ' + n ’ 


E vt — £ u a u “h £ v a v 


where, J? u and E v are constant. 


57 A 57 A 

V7 = u u + — a v 
5u at; 

Substituting (59) and (60) into the cost functional, (54), gives: 

NE * r BT BT 

F(T) = T' \( — d u + —d v ) + (E u d u + E v d v )\ 2 dudv 

e^i. Js* Ou Bv 

NE ‘ r BT BT 

= Ejj^+^ + ^ + E^aXdudv 

NE ‘ t BT BT 

= S X> l( lh + K)1 + { !h + E -^ iuiv 

The minimization of the cost functional is as follows: 

w- 

- Pl/I +<£> + <£- +<<£»** 

f ,,BT „.,dNi. ,dT _ w BN,.., , 

= E 2 X.K aU + + ®.K *7)]^" 

= 6 

where, i = m’, n’. 

Substituting (57) and (58) into (62) yields the following: 

BF(T) ™' n r 9N t dN m . BN. BN 

~b¥~= 

' , BN,> BN' _ diV< „ w c Wi„. . 

+ ( Bv Tl> + Bv Tm> + Bv n ' + Bu ^ U V 

r BNjdN? dNjdNf MBN m . B_NBN „ 

_ Js « Bu Bu Bv Bv 1 Bu Bu Bv Bv 

BN BN n > BNi BN ' . _ BN t _ BN, „ , 

^17 + *T^r )r » + !Z b ' + -fo E -V udv = 0 


n 



where i = l' ,m ,n 


In the previous equation, all partial derivative terms are constant because of the nature 
of the interpolation polynomials in first order elements, so that the elemental equation can 
be expressed on S' as follows: 


A(^ + ^)T, + A( 


dN,dN m . mdN_ ,■ 

I r\ a rr 


{ ~^-d^ + ~d^~d^ )T ^~ A{ ~d^ Eu+ dv E ' 


* • 

i = ( ) m,n 

where, A is the area of the triangular element, e. 


(64) 


The previous equation can be expressed in matrix form as follows: 


k'l 1 k'm 


fo 1 1 1 Ic 

m t 1 

\ k / i / fc ' / 

L n ( m 


tyn' 


> t k * t 

mm m n 



TV 1 





= 

Cj 


J 


LC'n' J 


where, 


kj = 


Ci = 
i,j = 


(WON, dNjdNj 

du pu dv dv 

\ fy dv 

l , m , n 


(65) 


( 66 ) 


4.2 Local Coordinate System 


The local coordinate system (d u ,a„,a„), see Fig.2.3, can be calculated as follows: 


„ . . . . 1 m X m'l' t . Vm'l 1 . , Z mY - 

a u = + Uydy + u 2 a z = -====; = t d x + 7 a y + 7 a 2 

lm| L, m 'i' L ml' 


(67) 


or, 


a u = 


u x 

Uy 

u z 


J m'V 


X m'l' 

Vm'l 1 

Z m'l' 


where 

’ X m'l' ‘ 

VmY 

= 

’ *m' 
Vm 

— X t ’ 

-Vl' 


. Z m'l' . 


. Z m' 

~ z l' . 


and L m ' t ' = yY 


( I'm ' x l'n') 1 

a w = w x a x + w y a y + w z d z = ; 7 r | — <) a 

|t m x / n | ■AAi 


£y 

+ <Y* 


a x 

dy 

h z 

X m't 

' Vth'I* 

Z m'i‘ 

X n'l' 

l tn’l* 

Z n l' 


( 68 ) 


( 69 ) 
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Wx 

1 

y m 'l' Z n'l' 

Vn l' Z m'l' 

a w — 

Wy 

_ 2A 

X n l' Z m'l' 

1 

H 

3^ 

Is* 


- W * - 

. X m’l'y n 'l' 



( 70 ) 


where, 


’ X n'l' 


x n X l' 

Vn'l' 

= 

y n ' - vi' 

. Z n'l' - 


1 

'c 

N 


and A is the area of the triangle I mn on surface S e 


2A = + (ZnY^mY “ ®m'/' Z n'/') 2 + (^rn'i'^n'/' “ *n'l' PmY ) 2 ( 71 ) 


a,, = 


Vx&X d” Vydy “f" V z dz ^ 

a x 


1 


2AL m '/' 


1/ /,/Z — 1/ /./Z /»/ X 'i'Z — X VZ , r / X tjtV 1 1* — X '} f V > i f 

y-m l n / ?n l m / n i m ( m l n l m l n i n t <7m l 


X 1 1> 

m J 


J/mV 


Z 


(72) 


or, 


<2i» — 




( X n'l' Z m 'l' - X m'l' Z n'l') Z m'l' ~ farnYVn'l' ~ X n' I'Vm' l')Pm' l' 
( X m'l' Vn'l' ~ X n l' Vra'i' ) X m'l' ~ (PmT Z nY “ Vn'l' Z m'l') Z m'l' 

L (Pm'i'Vl' - yn'l' Z m'l')ym'l' ~ “ X m’ l' Z n' l' ) X m ' l‘ 


(73) 


4.3 Transformation from Global to Local Coordinates 


In the global system, the surface elements are three dimensional but in the local system, they 
are two dimensional. We need a transformation from the global system to the local system. 
Assuming node l' is located on the origin point in the local coordinate, one can obtain the 
nodal transformation as follows: 


£ 


L . 



V x Vy v z 


1 

Xj — Xl t 


Vi - Vi 

- 

. Z t ” Zf 


t I f 


i = l } n 


(74) 


The tangential component of (V, x B g ) on the surface, S e , can be expressed as follows: 

E u = E vt ■ a u (75) 

E v = E vt ■ a v (76) 

In the local coordinate system, the shape functions become two-dimensional, hence [3]: 


1 


N i‘ = TKiPi' + <h‘ u + r /' v ) 

\ 

N m' = 7>^(Pm' + <lm' u + T m' V ) 

N n > = 7^(P„' + ?n' U + V U ) 


( 77 ) 
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and 


dNf 

a% 

du 


3L. 

2A’ 

9ml 

‘ 2A ’ at 


dNi> 

<fcj 


2A' 


JY_. 
2A’ 
Imi. 
' 2A ’ 

T’ 


dv 2A 


where 


ft' = ftn'ft»' “ U n' U m'; ft' = U m' ” «»' J ft' = U n' “ ftn' 

Pm’ = U n'ft' ~ ft'Vi ftn' = V “ ft' ! ftn' = ft' ~ U n' 
Pn' = ft'ftn' ~ U m' ft' > ft*' = ft' ~ ftn' 5 V = ftn' “ ft' 

Substituting (78) into (66), one obtains the following: 


(78) 


(79) 


K 

u 

where, 


C, 

ft; 

5 FE Formulation Using First-Order Tetrahedral El- 
ements 

5.1 Shape Function 



k,'> 

l m 

fc /V 

»'/' 

fc / j 

m m 

fc m'n' 

V 

fc f * 

m 

fc nV 



rr, 1 


Ct 


T 1 

m 

= 

Cm' 


L^n' J 


Cw \ 


(80) 


= 7^ (ft ft + ftft) 
= g.-Eu ft 

i' ' ' 

= / , m , n 


(81) 


Fig. 2. 4. shows a first order tetrahedral finite element by which we can construct an interpo- 
lation function as follows: 


f 4>(x, y , z) = ai + a 2 x + 0:31/ 4- a 4 z 
\ ft) = <f>i i = l,m,n,s 

On nodes, l,m,n,s , the interpolation functions are as follows: 


(82) 


ai 

+a 2 X[ 

+CX3VI 

+a 4 zi 

= <h 

« 1 

+a 2 x m 

-\-0-3ym 

+ <*4 ftn 

= 


+&2 x n 

+ct3 y n 

-\~cx 4 z n 

= 

<*1 

+OJ2X, 

+<*3y. 

+a 4 z. 

= 4>* 


( 83 ) 
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or, in matrix form: 



’ a 1 ' 


' 4>i 


OC2 


4>m 


Ot 3 


4>n 


. <*4 . 


. <t>s . 


1 xi Vi zi 

1 X m V m Z m 
1 X n 2/n Z n 

Ll®* y, z. 

The interpolation coefficients, ai, 02 , 0:3 and 04 , can be solved for as follows: 


o-i dm a n a a 

h b m b n b , 

C{ C m Cn C 3 

di d m d n d 3 


( 84 ) 


where, 


CK 1 ' 


1 Xl 

yi 

zi 
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Mi I 
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1 Xm 

Vm 

Zm 
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<2 3 


1 X n 

Vn 

Zn 


<t>rx 
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. 1 x. 

Vs 
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. <f>s . 


' a/ 

d>m 

o-n o a 


' 1 

Xl 

yi zi ] 

bi 

bm 

b n b a 


1 


Vm 


ci 

Cm 

Cn C, 


1 

X n 

Vn 

Z n 

. di 

dm 

I 

•> 

^3 


. 1 

x s 

Vs 

Z 3 . 



’ 4>i 


4>m 


4>n 


. 4>s . 


(85) 


Thus: 


<f)(x,y,z) = ai + a 2 x + a 3 y + a 4 z 
= {o>l4>l 4” ®m^m 4~ 

4'(£>/<^i + bmfim + b n <f> n + b,<f> s )x 
+ (ci<pl 4" Cmfim “1“ C n (^> n 4" O a (f> a )y 
+{di<f>i 4 " d m 4> m 4 " d n 4> n + d a <fi s ')z 


or. 


<j>(x,y,z) - Ni(f>i + N m (f> m + N n (f> n + N s <p, 
= (a/ 4- 4- c t y + diz)4>i 

4-(a m 4- b m x 4" CmJ/ 4" d m z^j<f> m 

4“(®n ”1“ b n x 4" C n y 4" d n z^<p n 
4 -(a, + b,x 4- c,y 4- d s z)<f> a 

Here, the shape functions are as follows: 

Ni = ( ai 4- bix 4- ciy + d t z) 

Nm — (n m 4" b m x 4“ Cm y 4~ d m z) 

N n = (a n 4- b n x 4 - c n y + d n z ) 

N s = (a, + b a x + c a y + d,z) 

where, 


( 86 ) 


(87) 


( 88 ) 
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1 ym Z m 


1 yi zi 

6, = - 

1 yn 

II 

$ 

1 y n z n 


l y, z a 


1 y, z. 


1 y/ 


l yi zi 

6n = - 

1 ym 
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1 ym z m 


l y* z* 
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X/ 
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1 X g Zg 
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1 Xj Zj 
1 X m ^m 
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/6V-, 
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1 x m J/m 

1 Xn J/n 
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1 X/ 

y/ 



1 Xm 

ym 

/6V 
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dN, 
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(89) 


(90) 
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5.2 FE Formulation 


Using the first order tetrahedral element, one can discretize elemental equation, (52), as 
follows: 



r dNidNj_ dNj dNj 

Jn> x dx dx dy dy 


+ cr z 


dNj dN, 

dz dz 


}<m 


= (cTxKb J + cryCiCj + a z d l d J ) • V 


k = 


N,ds + E 

J S f 


(93) 


&1 (V, X Bg) n — + 

1,J = l,m,n,s 


■ t jt ill 

j =1 ,m ,n , j 


S >T > 

*.j j 


t jt tit 

j ~L ,m ,n ,5 


where, A ; mn is the area of the triangle (nodes /, m,n); subscript, n denotes the normal direc- 
tion on surface SJ; symbol, j is the node on surface S { ; V is the volume of the tetrahedral 
element, e. 


6 Calculation of Current Density 


The current density is proportional to the electric field intensity. 
In region, 0^ : 


J\ — cr-iEi — o \ + E t ) — — <7i (V</>i — (U, x B g )) 


(94) 


where, a\ is the conductivity of region, Qj, which is homogeneous. 

From (47) and (48), the element potential and its gradient can be expressed as follows: 


*• = e N > m 

,=h ,= i ; 


(95) 


vr = E VJ V, - E VN,T, 


(96) 


j=k 
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The element current density can therefore be calculated by: 


A = -*i{(E - £ ViVjTj) - (7. x 5,)} (97) 

i=i[ 

where, ViVj (z = /,m,n,s) can be calculated using equations (89) through (92). 

In region, f&2 ; 

In this region, no (V x B) component and no potential jump are encountered. Hence, 
the current density of the element becomes: 

(98) 

i=l i 


7 Example 


In this work, an example (see Fig.2.5) is chosen to verify that the practicality and utility 

of the 3D-FE model developed in the previous sections. In Fig.2.5, we choose a cylinder as 

the moving object in place of the space station. Some geometric parameters and physical 

parameters are as follows: 

the length of the cylinder, L = 20 m 

the radius of the cylinder, R = 5 m 

the radius of the outer sphere boundary, Rq = 40 m 

the flux density along the negative z-axis, B g = 0.45G = 0.45 x 10 ~*Wb/m 2 
the speed of the cylinder along the negative y-axis, V, = 8km / sec. 
the conductivity inside the cylinder, <7\ = 5.7 x 10 7 1 / (f2 • m) 
the conductivity outside the cylinder, a-i = 5.7 x 10 5 l/(f) ■ m) 

The outer boundary condition in this example is assumed to be the natural boundary con- 
dition. 

The 3D-FE gird with a total number of elements, NE=3456, and a total number of 
nodes, NN=669, is shown in Fig. 2. 6. An FE grid with a total number of elements, NE=28128, 
and a total number of nodes, NN=4963, is chosen for the FE computation, whose surface 
elements on the cylinder and outer boundary are shown in Fig. 2. 7. The equipotential lines 
and the current density distribution along the plane of the cross section (z = 0, y > 0) were 
plotted from the results of this FE model see Fig.2.8 and Fig.2.9, respectively. The numerical 
results of the current density component J x vs radius on Plane x=0 are plotted in Fig. 2. 10. 
The numerical results of the scalar potentials located on the x-axis are plotted in Fig. 2.11. 
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8 Conclusions and Possible Future Work 


As a result of this investigation two technical papers on ballooning FE models and on twisted- 
pair transmission lines were presented at the IEEE INTERMAG-92 Conference, and were 
simultaneously published in the September issue on the IEEE Transactions on Magnetics 
(IEEE Trans, on Magn. Vol.30, No. 5, 1992). These papers are included here in Appendix 
(III) and Appendix (IV), respectively. 

Meanwhile, in order to obtain higher accuracy results regarding the present (V, x B g ) 
problem formulated in this report, the asymptotic boundary conditions (ABCs) should be 
adopted instead of the natural boundary conditions. High order ABCs can improve the 
accuracy of the numerical results. 

When the moving body is long and slender, a sphere shape outer boundary would 
be “storage memory intensive” and numerically inefficient. It would be highly desirable 
to choose an outer boundary that is conformal to the shape of the object. A formulation 
based on the ABCs should be derived for this 3D-FE scalar potential problem with an outer 
boundary of an arbitrary shape ( an ellipsoid is one possibility). 

In future work, the flux density distribution caused by the resulting pattern of induced 
current in the plasma and structure can be computed. This can be accomplished by a 3D- 
FE vector potential formulation using second order elements with the ABCs on the outer 
boundary, which should be designed to accommodate various shapes of the orbiting body 
under consideration. 
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Figure 2.1 Orbital motion of Space Station Freedom in the Earth’s geomagnetic 

field in the ionosphere. 
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Figure 2.3 The first order te 
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Figure 2.6 Geometry of the example. 
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Figure 2.6 Tetrahedral element grid for a cylindrical object with a sphere outer 
boundary, for NE = 3456, NN = 669. 
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Figure 2.7 FE surface grid for a cylinder object with a sphere outer boundary. 
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Figure 2.9 The current density distribution outside the object on the cross 
section (z = 0, x ^ 0, y ^ 0). 
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3 Radiation Sources 

This portion of the investigation centers on the development of computer 
models to calculate the radiation from simulated electromagnetic interference (EMI) 
sources on the space station. 

Section 3.1 briefly describes the plasma environment of the ionosphere within 
the altitude range corresponding to the orbit of the space station. It is shown that 
the anticipated frequencies of potential EMI sources fall well below the plasma 
frequency. Therefore, radiation at these frequencies will be attenuated by the 
plasma. Method-of-moments analyses of the radiation by sources on a straight wire 
and on a simple model of the space station demonstrate that both the currents in the 
structure and the electric fields around the structure may be calculated. Calculated 
results are also shown for the simple model of the space station embedded in a zero- 
order-plasma. This plasma approximation includes the charged particle density at 
the altitudes of interest, but neglects the geomagnetic field, thereby resulting in an 
isotropic rather than an anisotropic plasma. The results show that while the plasma 
attenuates the radiation, induced currents in the structure nevertheless produce 
radiated fields that exceed those that would be present in the absence of the 
conducting space station structure. Software that can be used to visualize contour 
diagrams of the calculated electric field distribution around the simple model of the 
space station is described. 

Section 3.2 describes the development of a wire-grid model of the space station 
that has been analyzed using the method-of-moments by the Electromagnetic 
Environmental Effects Laboratory at the Georgia Tech Research Institute. This is 
almost the same approach that the European Space Agency is using to predict EMI 
problems on their spacecraft. Calculations have been made for eight different cases 
comprising two source locations, two polarizations, and two frequencies, 20 KHz and 
5 MHz. The calculated results show top view and side view contour graphs of the 
electric field for each case - a total of sixteen graphs. Software that can be used to 
visualize both the wire-grid space station model and contour diagrams of the 
calculated electric field distribution around the space station is described. 

Section 3.3 briefly summarizes the available literature on plasma sheath 
waves. Excerpts from three recent publications, a journal article and two conference 
proceedings, are given. It has been shown (by others) that a region of low charge 
density, called a sheath, forms around a conducting object, like the space station, 
when it is embedded in the ionospheric plasma. A source on the conducting object 
can very efficiently transfer electromagnetic fields to distant locations by launching 
waves in the sheath region. The condition for transmission is the opposite of that for 
propagation in the ambient plasma as described in Section 3.1. Namely, sheath 
transmission is unattenuated for frequencies below the plasma frequency and is 
attenuated for frequencies above the plasma frequency. Sheath waves may therefore 
be an important factor in EMI analysis. 

Section 3.4 briefly summarizes the available literature on finite-difference time- 
domain analyses of electromagnetic wave propagation in anisotropic plasmas. The 
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abstracts of three publications, two journal articles and a conference proceeding, are 
given. These new analysis tools may be important for predicting the propagation of 
EMI signals in realistic models of the ionospheric plasma, i.e., a plasma model that 
includes both the charged particles and the geomagnetic field. 


3.1 Simplp Wire Model of Space Station Freedom in a Zero-order Plasma 


The space station will be embedded in the ionospheric plasma, a region of 
dense electric charge. As far as electromagnetic wave propagation in the ionosphere 
is concerned, the ionosphere can be modelled as an isotropic plasma if the 
geomagnetic field is neglected. Electromagnetic waves are then attenuated when 
their frequency is below the plasma frequency, but they propagate freely and 
unattenuated when their frequency is above the plasma frequency. Most importantly, 
the direction in which the wave is traveling is immaterial - exactly the same behavior 
occurs. However, when the geomagnetic field is considered, the plasma becomes 
anisotropic and the affect of the ionosphere on electromagnetic wave behavior 
depends on the direction of the electromagnetic wave relative to the direction of the 
magnetic field. All electromagnetic interference analyses to date have assumed an 
isotropic plasma, i.e., the geomagnetic field has been neglected. This assumption 
seems appropriate when what is needed is the basic behavior of electromagnetic 
interference in the presence of the plasma - to assume otherwise would greatly 
complicate any numerical calculations with little improvement in understanding the 
basic behavior of the interaction. 

The key parameter in studying electromagnetic wave propagation in a plasma 
is the plasma frequency 


a) 


p 


N 


Ne^ 


( 1 ) 


where N is the electron density in electrons per m 3 , e is the electron charge in 
Coulombs, ra e is the electron mass in kg, and eo is the permittivity of free space. 
Although the ionosphere consists of both electrons and ions, most of the interaction 
occurs with the electrons because they have a much smaller mass, and therefore are 
more responsive to the influence of electric and magnetic fields. That is the reason 
that the plasma frequency is written only in terms of the electron charge and mass. 

The variable factor in eq. (1) is the electron density N. Figure 3.1 shows the 
ionospheric charge density as a function of altitude for an equatorial orbit 
corresponding to the orbit of the space station. If only the electrons are considered, 
and an orbit altitude of 450 to 550 kilometers is assumed, it can be seen that log 10 iV 
falls in the range of 4.7 to 6.3 as shown in the outlined rectangular region in Fig. 3.1. 
This corresponds to an electron density N of 0.05 x 1CP to 2.0 x 10 12 electrons per m 3 . 

Figure 3.2 shows the calculated plasma frequency, f p = 2rcco p , for the assumed 
range of N. From the graph, it can be seen that the minimum plasma frequency of 
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2 Mhz corresponds to the minimum electron density. Since propagation occurs for 
frequencies above the plasma frequency and attenuation occurs for frequencies below 
the plasma frequency, attenuation is expected for all frequencies below 2 Mhz. A 
plane wave is attenuated by a factor of exp(-ar), where a is the attenuation factor in 
m 1 and r is the distance in m. The attenuation factor a is dependent upon the 
frequency of propagation co relative to the plasma frequency co p according to 


a 


2n 

TN 


0) 


\ 2 


V <*> 


- 1 


( 2 ) 


Figure 3.3 shows the calculated attenuation in dB/m as a function of frequency for 
different electron densities. The attenuation factor increases with increasing electron 
density. For a given electron density, the attenuation is fairly constant until the 
plasma frequency is approached, at which point the attenuation drops to zero 
according to eq. (2). 

The effect of the plasma on electromagnetic wave propagation is explicitly 
shown in Fig. 3.4. This shows the 0 - component of the electric field radiated by a 
short dipole as a function of distance away from the dipole for free space (zero 
electron density) and for electron densities from 0.01 x 10 12 to 1.0 x 10 12 per m 3 . It 
is clear that increased electron density results in increased attenuation. Here we 
have assumed that the near-field attenuation of the dipole fields occurs in a similar 
way as far-field attenuation. 

The above has considered the effects of the plasma on electromagnetic wave 
propagation, but without the presence of the space station structure. The influence 
of physical structure can be incorporated by using the method of moments, a well- 
known numerical electromagnetic analysis method. The method of moments can be 
used to find the currents and the local electromagnetic fields about a conducting 
structure which has been excited by a current or voltage source. The calculation 
finds both the currents and the fields. In fact, the solution procedure first finds the 
currents and then calculates the fields radiated by these currents. As an example, 
Figure 3.5 shows a 100 m wire (divided into 34 segments for analysis) excited by a 
short monopole antenna at the left end. The current distribution in the wire is then 
calculated for a 20 KHz, 100 KHz, 1 MHz, and 5 MHz source for the free-space (no 
plasma) case. The calculated current distribution is shown in Fig. 3.6. At the lower 
frequencies of 20 and 100 KHz, the normalized current distribution is the same. The 
smooth decrease of the current from the maximum at the source end to zero at the 
far end is characteristic of an antenna that is short relative to a wavelength. The 
wavelength at 100 KHz is 3000 m, so the structure is about 0.033 wavelengths in 
length at this frequency. The length of the wire in terms of wavelengths is even 
shorter at 20 KHz. As the frequency increases to 1 MHz, the distribution changes 
only slightly, as the structure is now about 0.3 wavelengths in length. However, at 
5 MHz, the structure is 1.67 wavelengths in length. Since we expect standing waves 
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to form with nulls every half wavelength, excitation at 6 Mhz should result in 
1.67/0.5 or 3.3 standing wave cycles along the wire length. Figure 3.6 therefore 
shows the expected result. There are also fields radiated by the wire, but they are 
not shown here. However, it is important to recognize that the currents in the wire 
and the fields around the wire are coupled and both will be excited by a source on 
this or any other conducting structure. 

Figure 3.7 shows a simple wire model of the space station for which near-field 
calculations have been made using the method of moments. The 100 m longtitudinal 
wire consists of 20 segments, with a one-volt source on the left-most segment. The 
source location was selected to simulate the switching associated with the 
photovoltaic solar panels. The two vertical wires, representing the habitation and 
laboratory modules, are each divided into 10 segments. The calculations account for 
the influence of the ionosphere by assuming a zero-order plasma environment. 
However, the effects of the plasma are incorporated only in the calculation of the 
radiated field and not in the calculation of the currents which produce the fields. 
Simple wire models like this are acceptable when the physical details of the actual 
structure are of a dimensional scale much less than a wavelength. At 1 MHz, the 
largest dimension of this model of the space station is only 0.33 wavelengths in 
length, so modelling with wires is justified. 

Calculations were first made for the model in free space. Figure 3.8 shows the 
results for a 1 MHz source located at the left end of the structure, corresponding to 
the origin (x = 0, y = 0) in Fig. 3.8. The contour lines are in units of dB pV/m. The 
results show the field strength to be relatively strong around the source as expected. 
The strong fields at the opposite end of the structure from the source indicate that 
conducted currents in the structure are an extremely efficient mechanism for 
transferring the influence of the source to other points on the structure. Based on 
the rapid decrease of the fields in the free-space region around the source, normal 
free- space considerations would predict fields at the far end of the stucture to be an 
order of magnitude less than what is calculated with the structure present. This is 
the reason why the conducting structure plays such an important role in 
electromagnetic interference analysis - it provides a current pathway for 
electromagnetic fields. 

The calculations were repeated using a zero-order plasma environment with 
an electron density N of 0.1 x 10 12 , corresponding to a plasma frequency of about 3 
MHz. These results are shown in Fig. 3.9 and the attenuating effect of the plasma 
is obvious as the fields around the space station structure are greatly reduced (recall 
that the units are dB pV/m). However, the fields very close to the structure are still 
significant. This is once again a result of the fields being transferred to other 
locations on the structure via the conducted currents. The method of moments 
solution for the plasma case did not use the plasma in the calculation of the currents 
(the structure current is the same in the calculation of Figs. 3.8 and 3.9), but only in 
the calculation of the fields radiated by the current. Were the plasma included in the 
current calculation, attenuation of the current would be expected and even lower 
fields close to the structure should result. Nevertheless, it is clear that the 
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conducting structure plays a significant role in producing fields close to the structure 
that are substantially greater than would be expected from consideration of radiation 
effects alone. 


Computer Visualization 

The calculated method-of-moments results in Figs. 3.8 and 3.9 may be viewed 
on a color monitor. An IBM-compatible PC with 612 Kb memory and (preferably) a 
hard disk is required. A graphics program called TOPO, licensed from Golden 
Software, Inc., may be used. 

First create a directory on the hard disk called TOPO and change to this 
directory. Then insert Disk #1 into the 3.6" floppy disk drive, switch to the floppy 
disk drive, change to directory TOPO, and copy all files to the hard disk (usually 
drive C:). Then switch back to the hard disk. The contour graphs may be generated 
as follows: 


type topo/cmd=fs (upper or lower case) 

when the program is loaded, press F2 to display Fig. 3.8 (free-space 
result) 

press ESC when finished viewing 

press ESC and RETURN to exit the graphics program 

repeat the above except type topo/cmd=zop to display Fig. 3.9 (zero-order 
plasma result) 


3.2 Method-of-Moments Wire-Grid Model of Space Station Freedom 

The method-of-moments may be used to find the conducted currents and fields 
radiated from sources on the space station. Early work used a simple wire model 
embedded in a zero-order plasma to determine the near-field electric fields. The 
model was then improved by more accurate modelling of the habitation/laboratory 
module and the addition of solar panels. The latest model is shown in Fig. 3.10. This 
model is constructed using a wire mesh with 720 wire segments. 

It is interesting to note that the European Space Agency is using the same 
electromagnetic analysis technique to develop EMC models of spacecraft (ESA report 
dated 8-28-91). Specifically, the method of moments is used to calculate the currents 
and the near electric field for a spacecraft excited by a voltage source. The method 
uses triangular surface patches rather than a wire grid to model the structure. The 
referenced report indicates that future work will include time-domain solutions, 
excitation by current sources as well as voltage sources, and the capability to include 
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dielectric materials as well as conducting surfaces. This work has apparently not 
included the effects of a surrounding plasma, probably because the applications 
currently being considered involve satellites in orbit above the ionosphere, in which 
case free space is a reasonable assumption. 

The wire-grid model of the space station has been analyzed by the 
Electromagnetic Environmental Effects Laboratory at the Georgia Institute Research 
Institute using the General Electromagnetic Model for the Analysis of Complex 
Systems (GEMACS) computer program, which is a commercial method-of-moments 
electromagnetic analysis program. The program assumes free-space conditions, so 
the presence of the ionospheric plasma is not accounted for. 

Four sets of calculations have been made, one for each of four different source 
configurations, two locations and two polarizations, as shown in Fig. 3.10 that are 
labelled #1, #2, #3, and #4. Two source frequencies, 20 KHz and 5 MHz, are 
considered for each source location and the distribution of the electric field is 
displayed in two planes, the x = 0 plane (Side view) and the z = 0 plane (Top view), 
as defined in Fig. 3.10. A total of 16 calculations have been made and graphs of the 
results are enclosed in this report as Figs. 3.11 - 3.26. The key to the graphs is 
shown in Table I. 

Each of the graphs shows contours of electric field strength in dB pV/m for a 
one Volt source at the locations shown in Fig. 3.10. The reference for the dB scale 
differs from that used in Figs. 3.8 and 3.9 in that the reference level (0 dB) is set at 
the source resulting in a relative dB scale, whereas the simple wire model results 
displayed in Figs. 3.8 and 3.9 were on an absolute basis. The only practical difference 
is that the contours displayed here show negative dB values. 

While calculations have been made for 20 KHz and 5 MHz (the 260th harmonic 
of 20 KHz), the 20 KHz results are of the most interest. The 20 KHz results for the 
source at location 1, top (file C20_T1) and side (file C20_S1) views, show the locations 
of high and low field strength. In particular, they show, as in the previous results 
for a simple wire model, that the presence of the conducting structure results in 
higher electric fields away from the source than would exist if the structure were not 
present. As previously explained, electric currents in the structure serve as a 
transfer mechanism for the electric fields. 

Computer Visualization 

The model shown in Fig. 3.10 may be viewed and rotated on a color monitor. 
An IBM-compatible PC with 512 Kb memory, (preferably) a hard disk, and a mouse 
is required. A freeware program called 3dv may be used. 

First create a directory on the hard disk called 3D and switch to this directory. 
Then insert Disk #1 into the 3.6" floppy disk drive, switch to directory 3D, and copy 
all files to the hard disk (usually drive C:). The model may be viewed as follows: 

type 3dv (upper or lower case) 
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use the mouse to click on the correct graphics mode (if necessary) 
use the mouse to click on the o bullet to the left of FILE 
use the mouse to click on SSF.3DV 

use the mouse to click on the o bullet to the left of SHOW 

rotate the model by moving the mouse 

when through viewing, click either mouse button 

use the mouse to click on the o bullet to the left of EXIT 


The 16 calculated method-of- moments results for the wire-grid model may be 
viewed on a color monitor. An IBM-compatible PC with 512 Kb memory, (preferably) 
a hard disk, and a mouse is required. A graphics program called GraphTools, 
licensed from 3-D Visions Corporation, may be used. 

First create a directory on the hard disk called GT and switch to this directory. 
Then create a subdirectory on the hard disk called TMP. Then insert disk #2 into the 
3.5” floppy disk drive, switch to the floppy disk drive, change to directory GT, and 
copy file GTOOL.EXE to the hard disk (usually drive C:). Then switch back to the 
hard disk and change to subdirectory TMP. Then switch back to the floppy disk drive 
and change to subdirectory TMP. Copy all 33 files in TMP to the hard disk (usually 
drive C:). Then switch back to the hard disk and change to directory GT. Run 
program GTOOL to create the GraphTools graphics package. When the program has 
finished running, delete GTOOL.EXE. The contour graphs may be viewed as follows: 

type gt (upper or lower case) 

from the menu, use the mouse to click on File 

use the mouse to click on Directory 

use the mouse to click on OK 

use the mouse to click on File 

use the mouse to click on Load 

press F3 to view the available files 

use the mouse to click on the file that is to be plotted (see Table I) 
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use the mouse to click on OK 
use the mouse to click on Redraw 

to view another file, click on File, Load, and press F3 as before 
use the mouse to click on Quit when finished 


3.3 Sheath Waves - A Summary of the Recent Literature 


J. J. Laurin, G. A. Morin, and K. G. Balmain, “Sheath wave propagation in a 
magnetoplasma,” Radio Science 24, pp. 289-300, 1989. (Complete paper in Appendix 
I.) 

“It has been understood by many authors that the presence of a low 
electron density sheath surrounding an antenna immersed in isotropic 
plasma provides a propagation region for electromagnetic waves, a 
situation that is of particular interest at low frequencies (signal 
frequency co < plasma frequency co p ) where the ambient plasma is cut off 
for uniform plane waves. In this case, the metal-sheath-plasma region 
becomes analogous to a waveguide in which surface waves, so-called 
sheath waves, may be excited and propagate over sufficiently long 
distances ... also provide a possible mechanism for locally generated 
electromagnetic waves to propagate between widely separated points on 
large spacecraft in the ionosphere.” 

G. B. Murphy and H. B. Garrett, “Interactions between the Space Station and the 
environment - a preliminary assessment of EMI,” Proceedings of the 3 rd Annual 
Workshop on Space Operations Automation and Robotics, Johnson Space Center, pp. 
493-507, March 1990. 

“If the structure is conductive it will have a significant sheath 
surrounding it due to the v x B motional potential. This sheath has 
been shown to be capable of conducting noise over large distances very 
efficiently. Sheath waves are guided waves that are conducted along 
conductors surrounded by sheaths much like waves in a coaxial cable 
transmission line. Anywhere sheaths overlay, the waves can propagate. 

The significance of this is that noise generated locally can be conducted 
along Space Station structures to other cables which may be sensitive 
to this frequency. Unless the sheath is forced to collapse, the waves 
propagate with little attenuation. Therefore, as a worst case scenario, 
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we assume that cables placed anywhere externally on the Space Station 
may be within a sheath which is connected to a source of noise via the 
structure-sheath coax transmission line.” 


K. G. Balmain, H. G. James, and C. C. Bantin, “Magnetoplasma sheath waves on a 
conducting tether in the ionosphere, with applications to EMI propagation on large 
space structures,” Proceedings of the 4 th Annual Workshop on Space Operations 
Automation and Robotics, Johnson Space Center, pp. 646-654, January 1991. 
(Complete paper in Appendix II.) 

“The evidence suggests that, on any large structure in low earth orbit, 
transient or continuous-wave electromagnetic interference, once 
generated, could propagate over the structure via sheath waves, 
producing unwanted signal levels much higher than in the absence of 
the ambient plasma medium.” 


3.4 Finite-Difference Time-Domain Modeling of the Ionosphere 

T. Kashiwa, N. Yoshida, and I. Fukai, “Transient analysis of a magnetized plasma 
in three-dimensional space,” IEEE Trans. Antennas Propag. 36, pp. 1096-1105, 1988. 

Abstract 

In recent years electromagnetic analyses of magnetized plasma have 
become important in various fields. Because of gyroelectric anisotropy 
with dispersion characteristics, the analysis is complicated, and three- 
dimensional treatment is indispensable. For efficient three-dimensional 
analyses in the time domain, the finite-difference time-domain method 
and the transmission-line-matrix method are generally used. However, 
until now these methods do not seem to have been applied to gyrotropic 
anisotropy with dispersion in the time domain. A new method was 
recently proposed for transient analysis in three-dimensional space 
based on the equivalent circuit of Maxwell’s equations and formulation 
by the Bergeron method. In this method, magnetized cold plasma, that 
is, gyroelectric anisotropy with dispersion, is formulated by adapting the 
trapezoidal rule to the characteristic equation of motion of the plasma 
in the time domain. The formulation of the magnetized cold plasma in 
three-dimensional space and in the time domain by the present method 
is described. The two principal problems, transverse and longitudinal 
propagation, are studied. By examining these results, the validity of 
this formulation is demonstrated. 
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R. J. Luebbers, F. Hunsberger, and K. S. Kunz, “A frequency-dependent finite- 
difference time-domain formulation for transient propagation in plasma,” IEEE 
Trans. Antennas Propag. 39, pp. 29-34, 1991. 

Abstract 

Computation of transient electromagnetic propagation through plasma 
has been a difficult problem. Frequency domain solutions are available 
for one or two dimensions and stratified plasma geometries. Transient 
solutions have been obtained via transformations to the tune domain, or 
by microscopic solution of the equations of motion of the charged 
particles. The finite-difference time-domain (FDTD) method is capable 
of explicitly computing macroscopic transient electromagnetic 
interactions with general three-dimensional geometries. However, 
previous FDTD formulations were not capable of analyzing plasmas for 
two reasons. First, FDTD requires that at each time step the 
permittivity and conductivity be specified as constants that do not 
depend on frequency, while even for the simplest plasmas these 
parameters vary with frequency. Second, the permittivity of a plasma 
can be negative, which can cause terms in FDTD expressions to become 
singular. A new FDTD formulation for frequency dependent materials 
(FD^TD has been developed, which removes the above limitations. In 
a previous paper (FD^D was applied to computation of transient 
propagation through a polar dielectric. In this paper, we show that 
(FD) 2 TD may also be applied to compute transient propagation in 
plasma when the plasma can be characterized by a complex frequency- 
dependent permittivity. While the computational example presented in 
this paper is for a pulse normally incident on an isotropic plasma slab, 
the (FDy^TD formulation is fully three-dimensional. It can accomodate 
arbitrary transient excitation, with the one limitation that the excitation 
pulse must have no zero frequency energy component. Time-varying 
electron densities and/or collision frequencies could also be included. 

The formulation presented here is for an isotropic plasma, but extension 
to anistropic plasma should be fairly straightforward. 

P. M. Franke and L. J. Nickisch, “Finite Difference - Time Domain solution of the 
Maxwell equations for the dispersive ionosphere,” Proceedings of the 1991 APS/URSI 
International Symposium, London, Ontario, p. 85, June 1991. 

Abstract 

The Finite Difference - Time Domain (FDTD) technique is a conceptually 
simple, yet powerful method for obtaining numerical solutions to 
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electromagnetic propagation problems. The application of FDTD 
methods to problems in ionospheric radiowave propagation is, however, 
complicated by the dispersive nature of the ionospheric plasma. In the 
time domain, the electric displacement is the convolution of the 
dielectric tensor with the electric field, and this convolution is 
untractable in the nearest neighbor FDTD approach. 

This difficulty can be avoided by returning to the constitutive relations 
from which the dielectric tensor is derived. By integrating these 
differential equations simultaneously with the Maxwell equations, 
dispersion (and absorption) is fully incorporated. The method is shown 
to be accurate by comparing to a known analytic solution for a special 
case. 
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Fig. 3.1 


Ionospheric charge density as a function of altitude. 





Fig. 3.3 


Attenuation vb frequency for specified electron densities. 
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Fig. 3.4 


Effect of varying the electron density on the electric field. 
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Fig. 3.6 Wire model with a short monopole antenna excited at one end. 


Ifonopole Source on End of 100 m Structure 



Fig. 3.6 


Current along a 100 m wire with a monopole source at one end. 



Fig. 3.7 


Simple wire model of Space Station Freedom. 







Fig. 3.10 (a) Wire-grid model of the space station, (b) Section of solar panel 

showing the location of the two sources and the two polarizations. 


Electric Field Intensity in Normalized dB 
at 5 MHz with Z = 0 (720 seg model/Source at Pos 1 ) 

(Max = 0.325 V/m) 
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Electric Field Intensity in Normalized dB 
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Electric Field Intensity in Normalized dB 



Pile C5M S2 



File 020 T2 



Electric Field Intensity in Normalized dB 
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File C5M_S3 


Electric Field Intensity in Normalized dB 
at 20 kHz with Z = 0 (720 seg model/Source at Pos 3) 

(Max = 0.00664 V/m) 
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Fig. 3.21 



Electric Field Intensity in Normalized dB 



Filp C9 , 0 S3 


Electric Field Intensity in Normalized dB 
at 5 MHz with Z = 0 (720 seg model/Source at Pos 4) 

(Max = 0.0817 V/m) 
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File C5M_S4 




File C20 T4 



Electric Field Intensity in Normalized dB 



Filp C20 S4 


Table 3.1 Kev to Wire-Grid Graphics Files 


File C5M_T1 
File C5M_S1 

6 MHz 
6 MHz 

Top view 
Side view 

Source at location 1 
Source at location 1 

File C20_T1 
File C20_S1 

20 KHz 
20 KHz 

Top view 
Side view 

Source at location 1 
Source at location 1 

File C5M_T2 
File C6M_S2 

6 MHz 
6 MHz 

Top view 
Side view 

Source at location 2 
Source at location 2 

File C20_T2 
File C20_S2 

20 KHz 
20 KHz 

Top view 
Side view 

Source at location 2 
Source at location 2 

File C6M_T3 
File C6M_S3 

6 MHz 
6 MHz 

Top view 
Side view 

Source at location 3 
Source at location 3 

File C20_T3 
File C20_S3 

20 KHz 
20 KHz 

Top view 
Side view 

Source at location 3 
Source at location 3 

File C6M_T4 
File C5M_S4 

6 MHz 
6 MHz 

Top view 
Side view 

Source at location 4 
Source at location 4 

File C20_T4 
File C20_S4 

20 KHz 
20 KHz 

Top view 
Side view 

Source at location 4 
Source at location 4 
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J. J. Laurin, G. A. Morin, and K. G. Balmain, “Sheath wave propagation 
in a ma gnetoplasma,” Radio Science 24, pp. 289-300, 1989. 
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A theory is presented for propagation of electromagnetic wave* concentrated is the vacuum 'gap 
sheath region separating a pUoe conductor from a cold mignetoplasma. The applied static magnetic 
field is paraJW to the conductor and ia the direction of wave propagation. The theory predicts 
propagation at all frequencies between zero and the upper-hybrid frequency divided by y/t Experi- 
mental results air reported in which resonances caused by the wave refiectioB at the cod of a thin 
cylindrical antenna aligned with the magnetic field are used to determine the sheath-wave dispersion 
relation. The effect of the antenna bias, plasma density and antenna diameter on the dispersion relation 
are studied experimentally and the results art compared with the plane-surface theory, la spite of the 
difference between the theoretical and experimental geometries, qualitative agreement is noted with 
respect to variations in magnetic field plasma density, sheath thickness and antenna radius. The 
prediction that the presence of an applied magnetic field raises the upper sheath wave cutoff frequency 
has been confirmed experimentally, as weD as the prediction that the presence of the sheath allows 
propagation above the cyclotron frequency. In particular, both theory and experiment show the exis- 
tence nor the cyclotron frequency of a special frequency at which the sheath-wave wave number is 
independent of the plasma density 


1 INTRODUCTION 

It has been understood by many authors that the 
presence of a low electron density sheath sur- 
rounding an antenna immersed in isotropic plasma 
provides a propagation region for electromagnetic 
waves, a situation that is of particular interest at low 
frequencies (signal frequency w < plasma frequency 
wj where the ambient plasma is cut off for uniform 
plane waves. In this case, the metal-sheath-plasma 
region becomes analogous to a waveguide in which 
surface waves, so-called sheath waves, may be excited 
and propagate over sufficiently long distances to sat- 
isfy resonance conditions set by the finite dimensions 
of the antenna structure. The immediate result of 
these resonances is a significant contribution to the 
antenna impedance. In fact, impedance measurement 
of a simple monopole antenna has probably been the 
most popular method used so fir in experimental 
sheath wave characterization [ Marec , 1970, Meyer ft 
al. 1974; hhizone et aL 1970c, b , c]. Sheath waves 
not only affect antenna performance but also provide 
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a possible mechanism for locally generated electro- 
magnetic waves to propagate between widely sepa- 
rated points on large spacecraft in the ionosphere. 
These waves could give rise to electromagnetic inter- 
ference if picked up by sensitive devices distant from 
the source. 

In this paper, results of theoretical and experi- 
mental research on sheath wave propagation in a 
magnetized plasma are presented, for the case of 
propagation in the direction of the static magnetic 
field B 0 . In earlier work, Mushiake [1964] derived a 
dispersion relation for surface waves along a thin, 
unsheathed wire parallel to B 0 . His equation predic- 
ted propagation cutoff for frequencies greater than 
the cylotron frequency a) t . Adochi [1977] proposed a 
transmission line model to calculate the impedance 
of a sheathless dipole antenna of finite thickness ori- 
ented at an arbitrary angle with respect to B 0 . In the 
case where the antenna and B 0 are aligned, this 
model predicts resonances called “lower n wave- 
length" below w t . Given the antenna length, the fre- 
quencies of these resonances can be predicted by Mu* 
shiake's formula. Miller [1967a, fc] found that the 
presence of a vacuum gap sbeath allows surface wave 
propagation 'therefore called sheath* wave propaga- 
tion) well above w f On the experimental side, 
measurements done in laboratory magnetized plas- 
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Fig. 1. Planar configuration for theatfe wave propagation in the /direction. 


mas [ Ishizone et aU 1969, 1970a, b ; Sawaya et al ^ 
1978] have revealed the presence of one, and some- 
times two, resonance peaks in the neighborhood of 
a> € . Due to this small number of resonances, there is 
only limited evidence that their presence is related to 
sheath-wave propagation. This is because there could 
be many other types of waves that arc likely to cause 
resonance in a finite-dimension plasma chamber. 
Also, even if the resonances reported were due to 
sheath waves, their small number would make diffi- 
cult if not impossible, the study of the dispersion 
relation. 

The objective of the present work is to establish 
clear experimental evidence that sheath waves indeed 
exist on an antenna in a magnetoplasma, that the 
waves can be easily excited, and that their dispersion 
relation can be determined approximately through 
measurement of antenna input impedance. A further 
objective is to relate the measured dispersion relation 
to the theoretical relation for sheath waves in a 
planar plasma-sheath-conductor region. 

I THEORY 

The derivation of the dispersion relation of surface 
waves propagating over an infinite met*] structure 
immersed in a magnetoplasma is now given. In order 
to simplify the problem from a mathematical and 
computational point of view, the following geometry 
has been adopted. The metal structure is reduced to 
an infinite plane over which is a semi-infinite mag- 
netoplasma A vacuum sheath of arbitrary thickness 
separates the metal from the plasma (see Figure 1). 

The study has also been limited to the particular 
case where the metal surface and the direction of 
propagation are parallel to the magnetic field which 


is in the i direction. In fact, it is the only situation in 
which the planar geometry is comparable to the cy- 
lindrical one in which we are actually interested In 
the general case, the angle between the normal to the 
surface and B 0 is not constant around a cylindrical 
surface, a situation that is impossible to reproduce 
with a planar structure. This incompatibility is 
avoided if the cylinder considered is aligned w ith B 0 . 

Finally, any variation of the field quantities in the 
y direction is neglected. This will not only simplify 
the analysis but also will provide similarity between 
the anticipated solution for the plane and the case of 
the cylindrical geometry in which azimuthal sym- 
metry would be assumed. 

For the plasma region, a cold magnetoplasma of 
uniform density was assumed. In this model, the rela- 
tive permittivity is given by the well known tensor; 



with 

X m io*/cj 2 Y « w r fw V * 1 — jv w 

w 9 electron plasma frequency; 
v) t electron cyclotron frequency 
o) excitation frequency; 

t effective electron-neutral collision frequency 


PRECEDING PAGE BLANK NO! FILLED 


LAURIN ET AL SHEATH WAVES IN MAGNETOPLASMA 


»| 


The SI system of units has been adopted in this 
paper. Also, all the field quantities, except B 0 , will be 
manipulated as phasors. 

In the case of harmonic time dependence of the 
type e hm , Maxwell’s equations are written 


V * H - ;'wt 0 K • E (la) 

V * E - H (2b) 

V • K • E - 0 (2e) 

V • H - 0 (2d) 

for the plasma region, and 

V * H ■ ;uk, E (3a) 

V * E - -ju>a 0 H (W 

V • E - 0 (3c) 

V • H - 0 (3d) 


in the vacuum-gap sheath region. Assuming a r vari- 
ation of the type t~** % where k is the complex propa- 
gation constant, one can write the V operator as 


in which fi\ « C, D, E and F are arbitrary 

constants. The boundary conditions at the metal sur- 
face (x = 0) require that £, * 0 and E, * 0. This sets 
£ * F « 0 in (4) and (5). 

2.2. Solution in the plasma 

Unlike the sheath region, the anisotropic nature of 
the plasma makes impossible the separation of the 
field components into two independent sets. There- 
fore it is only necessary to solve for one field compo- 
nent from which all the others can be derived. The 
choice of this field component may seem purely arbi- 
trary a priori. The choice of a component belonging 
to the TM mode will lead to a more general solution 
because, as it will be seen later, the components of 
the TE wave vanish if to, is reduced to 0. Also, in 
order to ensure conformity with the existing litera- 
ture, the solution for £, is found 

After a few algebraic manipulations [ Laurin, 
1986], equations (2) can be combined into the follow- 
ing: 

(Vr - - 0 <*) 


Oy. -ju) 

2.1. Solution in the sheath 

Using this definition of V, the components of E 
and H can be conveniently separated into two inde- 
pendent sets and the general solution to equations (3) 
is written as a superposition of a TM and a TE 
mode: 


TM mode 

H, - C cosh (Hx) + £ sinh (fix) 

(4a) 


k 

(4k) 


£ m — (C sinh {hx) + £ cosh (fcx)) 
* )ast 0 

(4c) 

TE mode 

E f m D sinh (Hx) + F cosh (Hx) 

(5a) 


-k 

H, E, 

<5fc) 


H m (D cosh (Hr) + F sinh (kx)) 

(5c) 

where 

h - (k 2 - fit)' 1 * 

(5/) 


where Vr is the component of V perpendicular to the 
i direction. The quantities aj and aj are solutions of 


the biquadratic equation, 

ox* + bs 3 + c * 0 pa) 

with 

e - K Pb) 

b - Pllli'* - X'* + XX 0 ) - k 2 (K 0 + XT (7c) 
c - foK 0 (K l - X* 2 ) - Ifilk’K,, K + k*K 0 Pd > 


With Vr * (d/cx)x. the general solution to (6) is 
E. - Ae‘" + Be ,, ‘ + Ce*"* + /»**'* 

A, B, C and / are constants. If we choose a, and i 2 
such that their real parts are negative, the solution 
for £. remains physically possibly only if C « / - 0. 
Thus, 

£, - Ae*'* + Be 1 " («) 

However, in the particular case where 

b* - *ae - 0 (9) 

the general solution for £.i s 

£. - [Ax * Bir" 0°) 

with 



292 


or this last solution, * can be obtained directly 
rom (9) and it is independent of the sheath thickness. 
However, calculation of the fields in the sheath and 
in the plasma using this value of k reveals that all the 
boundary conditions at the sheath edge cannot be 
satisfied simultaneously. This second type of solution 
should therefore be rejected. 

For the solution given by (8), the other field com- 
ponents are easily calculated by putting £ back into 
Maxwell’s equations. * ,nt0 
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(*%-») -g) 

€ — Bc“ iX 

a. 


e ~- J F!£*l 

H 

“ Aut, /C 0 e*’* +■ — e*” | 

(eMe*‘* + 0Be"*) 


hD cosh (hs) « (oAe*" + 0Be"*) (Jlrf) 

Note that in equations (lib) and (lid), the previous 
constant £> has been arbitrarily replaced by D'K’ 
This substitution facilitates the use of the same deri- 
vation for the solution in the isotropic case by simply 
setting a, * 0 in the expressions for the fields. When 

m " °’.r h * VC *' * K ° lnd ** « 0 and equations 
(7) can then be combined to form 


■] 


or 


which leads to 




9 m 8 m 0 


Also, the general solution for £, in the plasma region 
becomes 


H. - 


£, - At" 


kK’ wn 0 


( 12 ) 


where 


Equations (I la), (lib), and (lie) are then equivalent 
to, respectively, 


V « Kjk* -tlKl- Ki\ 
°=Ko(k'-e\K-)~K- 3 \ 

2.3. Dispersion relation for the sheath waves 
The dispersion relation is obtained from the 

Sr "° nd K i0nS ^ ‘ he Shcath -P^ma interface 
re ( ) absenoe of charge accumulation other 
than polarization charges and (2) absence of surface 
" rre "‘ (e ’ ec,nc and magnetic^ Those conditions are 
,s e ^ a * x “ 0) is continuous: 

~ cosh (hs) - JK t (J- r**' + i la) 

is continuous: 

D sinh (hs) - j ^ Ae , “ + — Be (i 
(3) E. is continuous : 


C cosh (hs) = tHhUl Af - 

a 

Dm 0 

C sinh (As) = ^2 Ae“ 


(l 3a) 

(13b) 
(1 3c) 


Elimination of the constants A and C gives the dis- 
persion relation for the isotropic case: 


Unb (hs) 


K„h 


lb) 


When there is no sheath, the boundary condition on 

mak« ,H "‘“I ?“*“ * ett A * 0 in (12 ^ w bich 
mak« all the fields vanish. This differs from the 

sheathless anisotropic case where the boundary con- 
dition on £, (given by (8)) sets A - - B, resulting in 
nonv&nishing fields, * 

For the nnisouopic cese, the eliminetion of the 
constants * ». C . end „ |e.d, , 0 ,h. dispenio. eel,! 

tIOD 


~~T c sinh (M • Ae“' + Be*** (]] C ) 

(4) H t is continuous (satisfied by case 2), (5) H is 
continuous (satisfied by case Ik and (6) H. is conti'nu- 


Unh‘ (hs)(9a t - + th ~ 


+ tinh (hs)fo 




(14) 


